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SUMMARY 

The viscous-shock-layer equations for hypersonic laminar nad 
turbulent flows of radiating or nonradiating gas mixtures in 
chemical equilibrium are presented for two-dimensional and 
axially-symmetric flow fields. Solutions are obtained using an 
implicit finite-difference scheme and results are presented for 
hypersonic flow over spherically-blunted cone configurations at 
freestream conditions representative of entry into the atmosphere 
of Venus. These data are compared with solutions obtained using 
other methods of analysis. 
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SYMBOLS 


damping factor [eqs. (15) and (16)] 


C. 

1 


mass fraction of species i, p^/p 


C Pri 


mass fraction of element % 


N 


frozen specific heat of. mixture, £ ^i C p i 


i=l 


specific heat of species i, C* ^/C* 


H 


H. 


U‘ 


defined quantity, « h + 

r 2 


total enthalpy, H + 
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N 


enthalpy of mixture, ^ C.h. 

i=.l 1 1 

enthalpy of undecomposed ablation material 


h. 

l 


enthalpy of species i, ht/U* 2 


flow index: 0 for plane flow; 1 for axisymmetric flow 


K 


thermal conductivity of mixture, K*/u* f C* 

r6i p f °° 


i 

M* 

M* 


mixing length [eq. (13) 3 


molecular weight 


molecular weight of mixture 


N 


number of species 


N 


Le 


Lewis number, p*D?.C*/K* 

xj P 
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turbulent Lewis number 


Prandtl number, y*C*/K* 

Jr 


turbulent Prandtl number, y*C*/K* 


Reynolds number, p*U*r*/y* 


Schmidt number, N gc = M p r / N Le 


coordinate measured normal to body, n*/r* 


normal coordinate [eq. (22)] 


pressure-gradient parameter [eq. (25)3 


pressure, p*/ p* 


( u ”) : 


divergence of the net radiant heat flux, Q*R*/p£u£ 3 

net radiant heat flux in n-direction, q*/p*U* 2 

' 00 00 

component of radiant flux toward the shock 
component of radiant flux toward the wall 
convective heat flux to the wall [eq. (13)] 


radius measured from axis of symmetry to point on body 
surface, r*/r* 


nose radius 


coordinate measured along body surface, s*/r* 


temperature, T */ T £ef 
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temperature, 2 /C* «> 

free-stream velocity 

velocity component tangent to body surface, u*/U* 
friction velocity [eq. (19)] 

velocity component normal to body surface, v*/U* 
scaled mean velocity component [eq. (18)], v /u, 

W L 

shock angle defined in figure 1 
angle defined in figure 1 
normal intermittency factor [eq. (30)] 
boundary-layer thickness 

incompressible displacement thickness [eq. (29)] 
number of atoms of the JLth element in species i 
normalized eddy viscosity, y T /]J 
eddy viscosity, inner law [eq. (12)] 
eddy viscosity, outer law [eq. (20)] 
transformed n-coordinate, n/n 

s 

body angle defined in figure 1 
body curvature 
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a 

a* 

< J > 1 / 2 1 3 

Superscripts : 

t 

j 0 for plane flow; 1 for axisymmetric flow 

quantity divided by its corresponding shock value 
* dimensional quantity 

' total differential or fluctuating component 

" shock-oriented velocity component (see fig. 1) 

Subscripts; 

e boundary- layer edge 

i ith species 


molecular viscosity# VVu*(T* ef ) 


eddy viscosity 


coordinate measured along body surface# E, s 


density of mixture, p*/p* 


Reynolds number parameter# 


»■* (T ^f) 

p*U*r* 

r oo oo fi 


1/2 


Stefan-Boltzmann constant 


quantities defined by equations (4b# 4c, 4d) 
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Jtth element 
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s shock 

w wall 

® free stream 

- values for the solid ablation material at the surface 

INTRODUCTION 

Numerical methods for calculating flow fields with ablation 
products injected into a radiating gas mixture in chemical equi- 
librium have been developed by Sutton (ref. 1) and Moss (ref. 2) 
for the analysis of hypersonic flow over blunt entry probe config- 
urations. The numerical solution procedure developed by Sutton 
is applicable to laminar, transitional, and turbulent flow and is 
obtained using a coupled inviscid flow-boundary-layer analysis. 

The inviscid flow solution is determined by an explicit time 
dependent finite-difference scheme similar to the method developed 
by Barnwell (ref. 3) , and the boundary-layer equations are solved 
by use of an integral matrix procedure (BLIMP) developed by Bartlett 
and Kendall (ref. 4) . 

Moss' analysis is restricted to laminar flow and solutions are 
determined using an implicit finite-difference scheme developed by 
Davis (ref. 5) for solving the viscous-shock-layer equations. The 
principal advantages of this technique are that the solution is 
direct and that the effects of inviscid-viscous interactions are 
included within a single set of governing equations which are 
uniformly valid throughout the shock layer. 

In the analyses presented by Sutton and Moss, the radiation 
heat transfer is calculated using the method developed by Nicolet 
(refs. 6 and 7) . This radiation model assumes a nongray gas and 
accounts for molecular band, atomic line, and continuum transitions. 

The present report presents the development of a viscous-shock- 
layer analysis applicable to laminar and turbulent flow of radiating 
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or nonradiating gas mixtures in chemical equilibrium. This analysis 
is based upon the viscous-shock-layer analysis applicable to tur- 
bulent flow of perfect gases developed by Anderson and Moss (ref. 8) 
and the laminar viscous-shock-layer analysis for equilibrium chem- 
istry developed by Moss (ref. 2) . 

Results obtained with the present method of analysis are 
compared with methods which include corrections for inviscid-viscous 
interactions. Solutions are presented for a 120-degree (total 
angle) spherically-blunted cone configuration at freestream condi- 
tions representative of entry into the atmosphere of Venus . 
Heating-rate distributions are compared for a cold wall (freestream 
temperature) nonradiating shock layer and a radiating shock layer 
with injected ablation products. 

The availability of comparative data obtained using methods 
corrected for inviscid-viscous interactions is limited. Consequently, 
the data obtained using the present method of analysis are to be 
considered as preliminary and serve primarily to establish stability 
of the numerical method. It is emphasized that no attempts have 
been made to obtain better agreement with either of the analyses 
used in the comparisons. A more extensive data base is necessary 
to establish the validity of the present solution procedure. 

ANALYSIS 

Governing Equations 

The equations of motion for reacting gas mixtures in chemical 
equilibrium are presented by Bird, Stewart, and Lightfoot (ref. 9) . 

The formulation of these equations in body-oriented coordinates 
appropriate for viscous-shock-layer analysis of laminar flow of 
radiating and nonradiating gases is presented by Moss (ref. 2) . 

For turbulent flow, the viscous-shock-layer equations are 
derived using methods analogous to those presented by Dorrance 
(ref. 10) for the turbulent-boundary-layer equations and are 
expressed in non-dimensional form for the coordinate system shown 
in figure 1 as: 
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Continuity 
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Energy: 
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Elemental continuity: 
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State: 


p = pTR*/M*C* w 
P f 


( 6 ) 
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Boundary conditions . The boundary conditions at the shock 

are calculated by using the Rankine-Hugoniot relations. At the 

wall, the no-slip and no-temperature- jump boundary conditions 

are used; consequently, u =0. The wall temperature and mass 

w 

injection rate are either specified or calculated. For the cal- 
culated mass injection conditions, the ablation process is assumed 
to be quasi-steady and the wall temperature is the sublimation 
temperature of the ablator surface. With these assumptions, the 
expression for the coupled mass injection rate is 


m 


-rr* 


— n * 


C , W 


‘r ,w 


"IT 

£ (C.h*) - 

1 1 w 

i=l 




(7) 


For ablation injection, the elemental concentrations at the wall 
are governed by convection and diffusion as given by the equation 



Precursor effects are neglected while the energy reradiated from 
the surface is included in the radiation transport calculations. The 
net radiative flux, q r , can be represented as the difference of 
two components 



(9) 


At the surface 


(+) * = ea*T* 4 


*r,w 


w 


where e is the emissivity of the ablator 


( 10 ) 
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The heat transferred to the wall due to conduction and 
diffusion is 


q 


c # w 




( 11 ) 


Radiative transport . The radiative flux, q^ f and the diver- 
gence of the radiative flux, Q, are calculated with the radiative 
transport code RAD, as presented in references 6 and 7. The RAD 
computer code has been incorporated in the present viscous-shock- 
layer computer code (HYVIS) 'and streamlined for computational 
efficiency . 

The RAD code accounts for the effects of nongray self-absorption 
and radiative cooling. Molecular band, continuum, and atomic line 

transitions are included. A detailed frequency dependence of the 

< 

absorption coefficients is used for integrating over the radiation 
frequency spectrum and the tangent slab approximation is used for 
integrating over physical space. 

Thermodynamic and transport properties . The equilibrium 
composition is determined by a free energy minimization calculation 
as developed in reference 11. Thermodynamic properties for specific 
heat, enthalpy, and free energy and transport properties for vis- 
cosity and thermal conductivity are required for each species 
considered. Values for the thermodynamic (refs. 12 and 13) and 
transport properties (ref. 14) are obtained by using polynomial 
curve fits. The mixture viscosity is obtained by using the semi- 
empirical formulae of Wilke (ref. 15) . 

Eddy-Viscosity Approximations 

A two-layer eddy-viscosity model consisting of an inner law 
based upon Prandtl's mixing-length concept and the Clauser-Klebanof f 
expression (based on refs. 16 and 17) for the outer law is used in 
the present investigation. This model, introduced by Cebeci (ref. 18) , 
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assumes that the inner law is applicable for the flow from the 
wall outward to the location where the eddy viscosity given by the 
inner law is equal to that of the outer law. The outer law is then 
assumed applicable for the remainder of the viscous layer. It is 
noted that the eddy viscosity degenerates to approximately zero in 
the inviscid portion of the shock layer. The degeneracy is expressed 
in terms of the normal intermittency factor given by Klebanoff 
(ref. 17) . The expressions used in the present investigation are 
given in the following sections. 

Inner-eddy-viscosity approximation . Prandtl's mixing- length 
concept is stated in non-dimensional variables as 


+ _ p£ 2 | 8U I 

e i a 2 u ' 3n ' 


( 12 ) 


The mixing length £ is evaluated by using Van Driest* s proposal 
(ref. 19) stated as 



where 



(14) 


Here, kj^ is the Von Kdrmdn constant, which is assumed to have a 
value of 0.4, and A + is a damping factor. 

Cebeci (ref. 18) suggests that for flows with a pressure gradient, 
the damping factor be expressed as 


+ + " */ 2 
A = 26(1 - 11. 8P) 


(15) 


and for flows with both a pressure gradient and mass injection, 

- 1/2 


\ + = 26 j- j^exp (11. 8v + ) - ij 


+ exp (11 


. 8v + )| 


(16) 
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where 


p + = - a 2 (* 2 -\ -4U 

X &) 0 

\ ' e t 


( 17 ) 


v + = ^ 

U T 


( 18 ) 


and 


u = o 


% / 3u\ 

P \ 3n / 


1/2 


( 19 ) 


Outer-eddy-viscosity approximation . For the outer region of 
the viscous layer the eddy viscosity is approximated by the Clauser- 
Klebanoff expression 


+ = k * pu e 6 lJi,n 
0 a 2 y 


( 20 ) 


where 


*■(* ( l ‘i) 


dn 


( 21 ) 


k, = O.Olt, 


and 


■i.n ‘ [ l * <f)‘] 


-1 


( 22 ) 
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Equation (22) is Cebeci's approximation (ref. 18) of the error- 
function definition presented by Klebanoff (ref. 17). 

For equilibrium flow without radiation, the boundary-layer 
thickness 6 is assumed to be the value of n at the point where 


H. 


H 


= 0.995 


t,°° 


(23) 


and is defined by linear interpolation in an array of local total 
enthalpies. This definition is approximately equivalent to the 
usual boundary-layer definition 



0.995 


(24) 


where u is the local value for the undisturbed inviscid flow 
e 

outside the boundary layer. 

The values of the parameters kj and k 2 in equations (13) 
and (20) depend on the flow conditions being considered, as does 
the constant represented by the value 26 in equations (15) and (16) . 
The values given are used for convenience in developing the numerical 
method. 


For radiating gases, the loss of energy from the shock layer 
makes the total enthalpy definition unsatisfactory. For these 
cases, the boundary-layer thickness, is assumed to be that portion 
of the shock layer which contributes 95% of the dissipated energy, 
and is defined by the expression 



V (1 + e ) 


dn 



y (1 + e ) 


0.95 


(25) 


dn 


This definition shows acceptable agreement with the total enthalpy 
definition when applied to nonradiating flows and should be an 
acceptable definition for radiating gases. 

Method of Solution 

Davis (ref. 5) presented a method for solving the viscous-shock- 
layer equations for stagnation and downstream flow. Moss (ref. 2) 
applied this method of solution to reacting multicomponent mixtures. 
The present method of solution is identical to that of references 
2 and 5. Therefore, only an overview of the solution procedure 
is presented here. 

The numerical computation is simplified by normalizing most 
of the variables with their local shock values. The transformed 
independent and dependent variables are 


g = n/n g 

p - P/P s 

K = s 

p = p/p s 

u = u/u s 

T = T/T 

s 

V = v/v s 

H = H/H 

s 

U = P/M s 

K = K/K 

s 




' (26) 
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Since the normal coordinate, n, is normalized with respect to 
the local shock stand-off distance, a constant number of finite- 
difference grid points between the body and shock are used. The 
transformations relating the differential quantities are 


n ’ 

Jl = A _£ A 

9s 9£ " n 0 3n 
s 


where 



dn 

s 


(27a) 


(27b) 


A = A A 

3n n g 9n 


(27c) 



and 


(27d) 


a 2 _ _i_ a 2 

9n ^ ” n 2 9n* 
s 

The three second-order partial differential equations are linearized 
and written in the standard form for a parabolic equation as 


9 2 W 

9n 2 


<*i 



a 2 W + 


a 3 


+ a 4 


9W 

3C 


0 


(28) 


where W represents tangential velocity for the s-momentum equation, 

K 

enthalpy for the energy equation, and elemental concentrations for 
the elemental continuity equations. For the energy equation, the 
divergence of the radiative flux is incorporated in the 0,3 term. 
When the derivatives in equation (28) are converted to finite- 
difference form by using Taylor's series expansions, the resulting 
equations are of the following form: 


AW n +BW + C W - = D 
n m,n-l n m,n n m,n+l n 


(29) 


The subscript n denotes the grid points along a line normal 
to the body surface, and the subscript m denotes the grid stations 
along the body surface. Equation (29) , along with the boundary 
conditions, constitute a system of the tridiagonal form and can 
be solved efficiently. 

A variable grid spacing is used in both the tangential and 

normal directions to the surface so that the grid spacing can be 

made small in the region of large gradients. The order of the 

truncation terms neglected are A£ (first order accurate) and 

either An An , or (An - An , ) . 

n n-1 n n-1 

The equations are solved at any body station m in the 
order shown in figure 2. The governing equations are uncoupled 
and the values of the dependent variables are computed one at a 
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time. Each of the second-order differential equations is individually 
integrated numerically by using the tridiagonal formalism [eq. (29)]. 
The global continuity equation is used to obtain both shock stand- 
off distance and the v components of velocity. By integrating 
the global continuity equation between the limits of 0 to 1, a 

quadratic equation for n is obtained. For the v component 

s 

of velocity at q, the global continuity equation is integrated 
with respect to q between the limits of 0 to q. The pressure, 
p, is determined at station m by integrating the normal momentum 
equation with respect to q between the limits of 1 to q. 

The equation of state is used to determine the density. The solu- 
tion is iterated until convergence is achieved. The solution 
advances to the next body station, m + 1, and uses the previous 
converged solution profiles as initial values for starting the 
solution at station m + 1. This procedure is repeated until a 
solution pass is obtained. 

The first solution pass provides a first approximation to the 

i 

flow field solution because the following assumptions are used. 

2 

The thin shock-layer form of the n-momentum equation, ■ - P— — • , 

' 1 + nK ' 

is used; the stagnation streamline solution is independent of 

downstream influence; the term dn /d£ is set to zero at each 

s 

body station; and the shock angle a is assumed to be the same 
as the body angle 0. These approximations are then removed by 
global iteration. Two solution passes are generally sufficient. 

This solution procedure is programmed for the CDC 6600 computer. 

DISCUSSION AND RESULTS 


Numerical solutions obtained with the present method of analysis 
are compared with an integral boundary-layer solution for a non- 
radiating shock layer with a cold wall boundary condition, and with 
a solution which couples the inviscid flow and boundary- layer 
equations for a radiating shock layer with surface ablation. 

For the nonradiating shock layer with a cold wall boundary 
condition, solutions were determined by Edquist* using the integral 

* The boundary-layer data have not been published and are presented 
by permission of C.T. Edquist, Martin Marietta Corp. , Denver 
Division, Denver, CO. 
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boundary-layer solution procedure (SHIV) discussed in reference 20. 
Freestream conditions correspond to a typical trajectory point for 
entry into the atmosphere of Venus. The freestream velocity, tempera- 
ture, and density are 10 km/s , 200 K, and 0.01 kg/m 3 , respectively. 

The atmospheric composition expressed in mole fractions is 0.95 C0 2 
and 0.05 N 2 . The body considered is a 120-degree (total angle) 
spherically-blunted cone having a nose radius of 0.368 m and a 
base radius of 0.66 m. The surface temperature is assumed to be 
equal to that of the freestream. 

The inviscid flow field solution used to specify edge conditions 
for the boundary- layer solutions was determined using a single strip 

integral method which accounts for the upstream influence of the 

< 

sonic corner. This influence cannot be accounted for in the 
present analysis, but as shown in figure 3, this influence is 
significant only in the region 1.6 _< r / r n £ 2. For r / r n < 1-6, 
the maximum difference in the surface-pressure distribution com- 
puted using the present method and the single strip integrated 

i 

method is less than 4%. 

Heat-transfer rate distributions corresponding to boundary- 
layer solutions for isentropic expansion edge conditions (and for 
edge conditions corrected for vorticity effects) are compared with 
the present method of solution in figure 4. It is noted that 
the present analysis was obtained assuming instantaneous transition 
from laminar to turbulent flow. The heat-transfer rate correlation 
formula used in the boundary-layer analysis includes a transition 
correction. Both the present analysis and the boundary-layer 
analysis corrected for vorticity effects show a significant 
increase in heat transfer when compared with the boundary-layer 
solution for isentropic expansion edge conditions. The present 
analysis and the corrected boundary-layer analysis differ by as 
much as 30% in the region of fully developed turbulent flow. 
Considering the assumption of local similarity used in the integral 
boundary-layer analysis and the different methods of turbulence 
modeling, the differences between the two methods of solution 
is not excessive. 
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The solution presented for a radiating shock layer with injection 
of ablation products corresponds to the Venus entry conditions used 
in the coupled inviscid flow-boundary-layer analysis presented by 
Sutton (ref. 1) . The atmospheric composition expressed in mole 
fractions is assumed to be 0.97 C0 2 and 0.03 N 2 . Freestream 
velocity, temperature, and density are 8.8 km/s, 180 K, and 0.0058 
kg/m 3 , respectively. The geometry considered is a 120-degree (total 
angle) spherically-blunted cone having a nose radius of 0.325 m and 
a base radius of 0.69 m. The ablator material is carbon-phenolic 
having a composition expressed in mass fractions of 0.11 O, 0.004 N, 
0.851 C, and 0.035 H. 

The surface ablation-rate distribution used in the present 

< 

analysis corresponds to that determined by Sutton (ref. 1) , and 
is shown in figure 5. Comparisons of the surface-pressure distri- 
butions and shock shapes corresponding to the present analysis and 
that of reference 1 are shown in figures 6 and 7, respectively. 

The maximum differences in the surface-pressure distributions and 
shock-layer thicknesses determined by the two methods of analysis 
are approximately 3%. Since the inviscid solution does not account 
for displacement effects, the differences noted are expected for 
the specified injection rates. Other properties within the 
essentially inviscid portion of the shock layer show similar 
agreement. 

Comparisons of radiative and convective heating-rate distri- 
butions corresponding to the two methods of analysis are shown 
in figure 8. Differences of 5 to 10% are obtained for the radia- 
tive heating-rate distributions, and convective heating-rate 
distributions differ by 10 to 15% in the region of laminar flow. 

In the turbulent flow region, the agreement between the two methods 
of analysis is unsatisfactory. The reason for the opposite trends 
in the turbulent heating-rate distributions corresponding to the 
two methods of analysis has not been determined. Additional calcu- 
lations will be necessary to resolve these differences. 
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CONCLUDING REMARKS 

The results of the present investigation demonstrate that 
humerically stable solutions to the viscous-shock-layer equations 
can be obtained for turbulent flows of radiating and nonradiating 
gas mixtures in chemical equilibrium. Acceptable agreement between 
the present method of analysis and an integral boundary- layer 
analysis is obtained for a nonradiating shock layer without injec- 
tion of ablation products. The agreement between the present 
method of analysis and a solution which couples the inviscid 
flow and boundary- layer equations is unsatisfactory for the case 
of a radiating shock layer with ablation products injected into 
the layer. The limited availability of comparative data obtained 
with other methods of analysis is not sufficient to verify the 
present method of analysis. 
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Figure 1. - Coordinate system. 
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Figure 2. - Flow design of local solution procedure 
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Figure 4. - Comparisons of heating-rate distribution 
on 120-degree sphere-cone. 
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Figure 5. - Carbon-phenolic ablation-rate distribution for 120-degree sphere-cone. 
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Figure 6. - Comparison of surface pressure distribution 
on 120-degree sphere-cone. 
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Figure 7. - Comparison of shock-layer thickness distribution 
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